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ABSTRACT 

We show that the source population of long gamma-ray bursts (GRBs) has declined by at least a factor 
of 12 (at the 90% confidence level) since the early stages of the Universe (z ~ 2 — 3). This result has 
been obtained using the combined BATSE and Ulysses GRB brightness distribution and the detection 
of four GRBs with known redshifts brighter than 10 52 erg s _1 in the 50 - 300 keV range at their peak. 
The data indicate that the decline of the GRB source population is as fast as, or even faster than, the 
measured decline of the star formation rate. Models for the evolution of neutron star binaries predict a 
significantly larger number of apparently bright GRBs than observed. Thus our results give independent 
support to the hypernova model, which naturally explains the fast decline in the progenitor population. 

Subject headings: gamma-rays: bursts - methods: data analysis 



1. INTRODUCTION 

The cosmological evolution of GRB progenitors at red- 
shifts z < 2 can, in principle, reveal their nature. Indeed, 
we have unambiguous star formation data (hereafter SF; 
see Porciani & Madau 2001 and references therein) for the 
declining stage which started after z ~ 2, which we can use 
as a reference evolutionary curve. If GRB progenitors fol- 
low this curve or decline even faster than it, then we have 
to conclude that GRBs are most probably associated with 
the collapse of supermassive stars (hypernovae, or "failed 
supernovae" as originally suggested by Woosley 1993; see 
also Paczynski 1998 and MacFadyen & Woosley 1999). If 
the decline of GRBs is slower than the SF decrease then the 
coalescing neutron star binary model would be supported, 
as it naturally provides a delay between star formation and 
bursts. 

GRB afterglow observations provide three lines of evi- 
dence in favor of the hypernova model as an explaination 
for the long duration GRBs. First, a large fraction of 
the afterglows are found near the central regions of their 
host galaxies (Bloom, Kulkarni & Djorgovski 2001). Sec- 
ond, features have been found in the light curves of three 
afterglows which can be interpreted as a supernova com- 
ponent (e.g. Bloom et al. 1999, Lazzati et al. 2001). And 
third, absorption features in some x-ray afterglow spectra 
(Galama & Wijers, 2001) and an emission iron Ka line 
(Piro et al., 2000) indicate a high metal column density 
along the line of sight. For more details see the review 
of Meszaros (2001) and references therein. Although none 
of these facts constitutes a decisive argument by itself, to- 
gether they strongly favor the hypernova model. However, 
the cosmological evolution of GRB sources can provide a 
new and completely independent test for the nature of 
GRB progenitors. 

The problem of deriving the GRB source evolution from 
the data is not simple and cannot be solved by a straight- 



forvard cosmological fit to the log N - log P distribution 
with an unknown GRB luminosity function. Despite the 
wealth of statistics on GRBs accumulated by the Burst 
and Transient Source Experiment aboard the Compton 
Gamma-Ray Observatory (BATSE) (see Fishman et al. 
1989), the bright end of the distribution still contains too 
few events to provide a conclusive \ 2 fh. For a review of 
cosmological fits to the log N - log P distribution see, e.g., 
Bulik (1999). 

Stern, Tikhomirova & Svensson (2001 - hereafter STS), 
using the redshift data and the BATSE GRB sample, 
demonstrated the cosmological decline of GRB progeni- 
tors. A fast decline, similar to that of SF or even faster, 
was preferred by the data, but the statistics were insuffi- 
cient to distinguish between the predictions of the hyper- 
nova and neutron star (NS) binary models at a significant 
level. 

In this work we incorporate the data of the Ulysses GRB 
experiment, which more than double the number of strong 
GRBs and allow a more reliable interpretation of the red- 
shift data. The main objective is to achieve a scientifically 
meaningful constraint on the NS binary model. 

In §2 we describe the data and the procedure used to 
cross- calibrate the Ulysses and BATSE GRB samples. In 
§3 we outline the model fitting procedure, including the 
cosmological model, the parameterized luminosity func- 
tion, and the GRB source evolution hypothesis. In §4 we 
present our results and show that the data require a very 
fast GRB progenitor decline, seriously challenging existing 
NS binary models. 

2. THE DATA 



We have used three independent data sets. The first 
contains 3255 BATSE GRBs with durations longer than 1 
s, found by Stern et al. (2000, 2001) in the off-line scan 
of the 1.024 s time resolution BATSE continuous daily 
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records for the entire 9.1 yr BATSE mission 6 . This is the 
largest essentially uniform GRB sample, and its efficiency 
matrix has been measured using a test burst method. The 
second data set is the Ulysses sample, consisting of only 
bright GRBs, which are the most important ones for the 
aim of the present work. The Ulysses GRB detector has 
amassed well over 10 years of data to date, and since the 
detector is in interplanetary space and is neither Earth- 
nor spacecraft-occulted, it has « 4ir sr sky exposure and 
a larger effective duty cycle than BATSE (useful data are 
recovered for more than 95% of the mission), thus more 
then doubling the number of bright GRBs. The Ulysses 
GRB data on over 800 bursts have appeared in eight cata- 
logs so far (Hurley et al. 1999a,b; Laros et al. 1997, 1998; 
Hurley et al., 2000a, b,c; Hurley et al. 2001a); the instru- 
ment description may be found in Hurley et al. (1992). 
The third data set consists of the GRB redshift data, or 
more specifically the data on the four intrinsically bright- 
est events out of 23 GRBs with measured redshifts (up to 
November 2001) 7 . 

The first two data sets were combined to form a single 
log N - log P distribution, i.e. the number of events versus 
the apparent peak brightness, P a , while the third data set 
was used to constrain the bright end of the hypothetical 
intrinsic peak brightness (Pi) distribution (the luminosity 
function). 



rate versus BATSE flux is shown in figure 1. 
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2.1. The BATSE and Ulysses samples, their 
cross- calibration, and the joint log N - log P 
distribution 

The BATSE sample includes 3255 "long" (duration > 1 
s) GRBs selected from 3906 GRBs in the sample of Stern 
et al. (2001), with the requirement that the counts in the 
second highest time bin exceed 50% of the counts in the 
highest time bin. This sample corresponds to 9.1 years of 
GRO observations with an average exposure factor 0.47. 
The latter includes both Earth-blocking and the average 
duty cycle for useful 1.024 s continuous records. 

The Ulysses sample covers 10.3 years (Dec 1990 - Feb 
2001) which overlap the entire BATSE mission. The aver- 
age exposure factor of the GRB experiment is ~ 0.96. This 
includes data outages, as well as periods when the back- 
ground was high due to solar proton events. The short 
events were removed from the sample using the same cri- 
terion that we applied to BATSE GRBs. The effective 
energy range of the Ulysses GRB detectors is w 25 - 150 
keV, while the BATSE data used for this study are for the 
50 - 300 keV band. 

Although both the Ulysses GRB detector and BATSE 
have quasi-isotropic angular responses, a direct conversion 
of Ulysses counts to BATSE peak flux is not possible for 
several reasons. First, the responses of both experiments 
have a weak angular dependence. Second, they operate in 
different spectral bands. And third, they have different ef- 
ficiencies as a function of energy. The only realistic way to 
construct a joint log N - log P distribution is to do a cross- 
calibration using joint Ulysses /BATSE events. There are 
278 such GRBs down to the Ulysses cutoff adopted here of 
100 counts s" 1 above background (the average background 
rate is 480 counts s _1 ). A scatter plot of the Ulysses count 



Fig. 1. — Ulysses count rate versus BATSE peak flux for 278 
GRBs detected by both experiments. The sample has been trun- 
cated at a Ulysses peak count rate of 100 s _1 . The dashed line 
shows the ratio of the Ulysses count rate to BATSE flux, 24.6 cm 2 , 
which we use as the calibration coefficient. 

The tendency towards higher Ulysses/BATSE ratios in 
the range below 10 photons cm~ 2 s _1 is caused mainly by 
the Poisson bias (i.e., when one selects the highest Pois- 
son fluctuation as a peak) and possibly by the hardness- 
brightness correlation in GRBs (see below). At the bright 
end of the distribution a similar bias is evident: the Ulysses 
count rates are again systematically higher. This is prob- 
ably a saturation effect in the BATSE count rates, result- 
ing from slow light emission in Nal scintillator (Meegan 
& Preece, 2001). This is probably a saturation effect in 
the BATSE count rates caused by dead time in the large 
BATSE detectors. 

To avoid these biases we have restricted the analysis 
to bursts in the 10 - 40 photons cm~ 2 s^ 1 range for the 
Ulysses/BATSE calibration. There are 73 events in this 
range; the average Ulysses/BATSE ratio is 24.6 cm 2 , and 
the rms variance is 5.2. The relatively large variance is 
due to the three reasons cited above. Note that this flux 
range gives the smallest Ulysses /BATSE ratio and there- 
fore the largest values of Ulysses peak fluxes. Therefore 
it provides a conservative constraint on the decline of the 
GRB population. 

In principle there should also be a brightness de- 
pendence in the Ulysses /BATSE ratio caused by the 
brightness-hardness correlation in GRBs (Nemiroff et al., 
1994) and the different spectral bands and efficiencies as a 
function of energy. This effect may contribute to the ten- 
dency towards higher Ulysses /BATSE ratios for weaker 
GRBs (figure 1). We cannot separate this from the Pois- 
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son bias. However, Atteia (2001) parameterized the cor- 
relation between the Ulysses count rate and the 50 - 300 
keV photon flux as as P a oc C 1 , . Thus the brightness 
dependence is weak and has a negligible effect on the re- 
sults. We prefer to use a constant calibration coefficient 
because there are insufficient data to quantify this depen- 
dence more accurately. 

The joint BATSE- Ulysses log N - log P distribution is 
shown in figure 2. It includes 77 Ulysses bursts above 
log(P a ) = 1.2. We chose this relatively high cutoff to 
conservatively avoid the Poisson bias in the Ulysses peak 
flux estimate. The number of BATSE GRBs in the range 
log(P ) > 1.2 is 43 (the total number of BATSE GRBs in 
the distribution is 3255). 
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Fig. 2.— The joint BATSE - Ulysses log N - log P distribution. 

The distribution below log(P a ) = 1.2 is normalised to 
the total BATSE exposure, 9.1 yr xF b , where F b = 0.47 is 
the BATSE exposure factor. The normalisation of points 
above log(P Q ) = 1.2 is the sum of the total BATSE/ Ulysses 
exposure during 9.1 yr and the Ulysses -only exposure dur- 
ing 1.2 yr: 9.1yrx(P b + P„-P b -P tl )+ 1.2yrxF u , where 
F u = 0.96 is the Ulysses exposure factor. 

2.2. The rate of intrinsically strong GRBs 

The sample of events with known redshift is subject to 
strong selection biases and cannot be used directly to de- 
termine the luminosity function (see STS). It does, how- 
ever, give useful information about the existence of intrin- 
sically very bright GRBs. We can use this fact to constrain 
the bright end of the hypothetical luminosity function: the 
predicted rate of GRBs with P, above some threshold at 
all redshifts should correspond to the observed rate. This 
constraint will affect the predicted number of apparently 
bright GRBs and therefore constrain the GRB source evo- 
lution model. 



For convenience we measure the intrinsic peak bright- 
ness Pi as the peak flux normalized to z = 1, taking the 
K-correction into account. The approximate relation be- 
tween Pi and the isotropic peak 50 - 300 keV luminosity 
is L = 3 • 10 50 erg s _1 P a . STS have chosen the range 
Pi > 40 photons s _1 cm~ 2 as a reference for the com- 
parison between the hypothetical luminosity function and 
the data, and we use the same threshold here. Four of 
the intrinsically brigthest bursts are above this threshold: 
GRB990123, GRB991216, GRB000131 and GRB010222 
, with redshifts of 1.6 (Djorgovski et al. 1999), 1.02 
(Vreeswijk et al. 1999), 4.5 (Andersen et al. 2000), 
and 1.477 (Stanek et al. 2001) respectively; their appar- 
ent fluxes, P a , in photons s _1 cm -2 , are 16.4, 67.5, 6.3 
(BATSE catalog) and 22.4 (Ulysses data with our Ulysses 
/BATSE calibration) respectively Their P, values are 45, 
69, 84 and 48 photons s _1 cm~ 2 , respectively. 

These 4 events were detected over 4.2 years from the 
beginning of 1997 to March 2001. The latter date corre- 
sponds to the end of the processed Ulysses data that we 
use here (currently there are three GRBs with measured 
redshifts after March 2001). 

The corresponding rate of GRBs at Pi > 40 photons s _1 
cm -2 , /40, multiplied by the probability that the burst will 
be detected and localized, its afterglow observed and its 
redshift measured (hereafter the sampling factor, F s ) is 



0.95. 



-0.60 
-0.38- 



Now, in order to estimate I40, we have to eval- 



uate the sampling factor. 

The first approach is simply to estimate F s for appar- 
ently bright events, which would give a reasonable up- 
per limit to the sampling factor for intrinsically bright 
events. Taking all Ulysses GRBs with peak count rates 
above 370 s^ 1 which corresponds to approximately 15 pho- 
tons s _1 cm -2 (62 events from 1997 January 1 to 2001 
March 1) we find redshift data for 4 of them (two of which 
are in the above list of four intrinsically bright GRBs). 
Using these numbers we estimate the sampling factor as 



0.064 



+0.038 
-0.026- 



Taking the ler upper limit, 0.1, as a con- 



servative estimate we obtain the rate of intrinsically bright 
GRBs I m ~ 10. 

The alernative approach is a direct estimate of the effi- 
ciency of the detectors used for burst localisations. 13 out 
of 23 redshift measurements were done using Beppo-SAX 
localizations. 7 of the remaining 10 bursts were localized 
by the Interplanetary Network (IPN: Hurley et al. 2001b). 
It is difficult to estimate the IPN efficiency, but the Bep- 
poSAX efficiency has a well defined upper limit: the total 
ficld-of-view of the two Wide Field Cameras, ~ 0.08 of 
the sky 8 . Not every localisation, even of a strong burst, is 
followed by the observation of an afterglow and a redshift 
measurenent. With this upper limit on the BeppoSAX 
efficiency and its share in the redshift sample we again 
obtain an estimate of the sampling factor of ~ 0.1. 

Therefore we adopt the estimate I40 = 10/yr as our 
baseline and, to take the poor statistics into account, we 
also rederive all our results for /40 = 4. Future observa- 
tions will show which value is closer to reality. 

3. FITTING MODELS 

The fitting model consists of three independent compo- 
nents: the cosmology, the evolution of the GRB source 
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population, and the intrinsic luminosity function (here- 
after just luminosity function or LF). 

The cosmological model is not very important for the 
purpose of the present work as it affects only large red- 
shifts while the main issue we are concerned with here 
is the source evolution at low redshifts. We adopted a 
flat vacuum-dominated cosmology (Oa = 0.7, Am = 0.3) ^ 
which is supported by recent data (see, e.g., Lukash, 2000). H 

The evolution of the source population is the objective § 
of our study. We tested four evolutionary functions. The ^ 
first is a non-evolving population (NE). The second is the t 
star formation function, which is a reasonable hypothesis 5 
for the evolution of GRB progenitors if they are collapsars. a, 
Porciani & Madau (2001) suggest three parameterized ver- ^ 
sions of the star formation rate which are very similar at ^ 
small redshifts, but differ at large redshifts where the in- 
terpretation of the data is ambiguous due to the poorly \ 
known effects of dust absorption. Again, the evolution at :§ 
large redshifts is beyond the scope of this work and we 
considered only one of these versions, namely: 



R Sf{z) = ,3.4,+' 22) Moyr^Mpc- 3 (1) 



This expression describes a constant SF rate at large 
redshifts and corresponds to a comoving volume. 

The two other evolution functions used correspond to 
neutron star merger models. We obtained them by con- 
volving the above SF rate with two different distributions 
for the delay between the formation of a binary system 
and the coalescence of its daughter neutron star binary. 
The first delay distribution was taken from Lipunov et 
al. (1995), hearafter L95, and the second from Portegies- 
Zwart & Yungelson (1998), hereafter PZY98. These distri- 
butions are quite different from one another. L95 predicts 
a peak at delays of 10 - 20 Myr and a long tail with a 
comparatively high probability of several Gyr delays. The 
distribution of PZY98 has a maximum around 1 Gyr and 
lower probability at several Gyr. 

The standard candle log N - log P distributions for these 
four models are shown in Figure 3. 

In addition to four fixed evolution models we tested dif- 
ferent slopes of the decline phase of the source population, 
modifying Eq. (1) as 



e 1.086in(a+l)2 
RSF ^ K ( e 1.086in(a+l)s + a ) ^ 



where a is a parameter describing the fall-off with red- 
shift: (a+1) is the ratio of the comoving rates of GRB 
emission at large z and at z = 0. The expression coincides 
with (1) at a = 22. 
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Fig. 3. — Standard candle log N - log P distributions for differ- 
ent models. The standard candle brightness corresponds to a peak 
flux of 1 photon s~ 1 cm' 2 at z = 1. Solid curve: SF model; dotted 
curve: NE model; dashed curve: SF model convolved with L95 delay 
function; dash-dotted curve: SF model convolved with PZY98 delay 
function. The crosses represent the observed log N - log P distribu- 
tion of 3255 long BATSE GRBs. Theoretical curves are normalized 
to the same integral number of GRBs while data points are given in 
an arbitrary normalization. 

The third component of the model is the hypothetical 
luminosity function. The data allow a wide choice with 
only two constraints: the width of the function, which 
must be at least 2.5 orders of magnitude (the luminosity 
range of GRBs with measured z), and the number of in- 
trinsically bright GRBs (see section 2) . We chose a broken 
power law, which proves sufficient freedom with a reason- 
ably small number of free parameters: 

dN/dP = C ■ P a - X for P 1 < P < P b , dN/dP = 
Ci • P^+P- 1 for P b < P < P 2 and dN/dP = beyond 
the interval [Pi, P 2 ]. The free parameters are a, (3, Pi, Pb, 
and C, while P 2 was fixed to a value above the maximum 
observed GRB brightness. 

We used the forward folding method when fitting GRB 
data, i.e., the hypothetical brightness distribution was con- 
volved with the efficiency matrix (1) and fitted to the ob- 
served differential log N - log P distribution (crosses in 
Figure 2) represented by 29 data points below P = 50 pho- 
tons s _1 cm~ 2 In 9.1 years of BATSE and Ulysses data, 
there were 15 GRBs brighter than this. We treat the range 
P > 50 photons s _1 cm~ 2 separately, estimating the like- 
lihood function of the fit for each peak flux range. For the 
main interval, this is the standard \ 2 probability function. 
For the tail of the brightness distribution, the likelihood is 
the Poisson probability of finding no more than 15 events 
apparently brighter than 50 photons s _1 cm~ 2 assuming 
an average number A§o predicted by the model. The final 
likelihood finction is the product of these two factors. 
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4. RESULTS 

The best fit integral log N - log P distributions for the 
four models, SF, PZY98, L95, and NE are shown in figures 
4 and 5. 
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Fig. 4. — Comparison between the predictions of different evo- 
lutionary models and the data. Solid curve: SF; dotted curve: NE; 
dashed curve: L95; long dashed curve: PZY98; straight dashed line: 
the Euclidean -3/2 slope. Crosses are the observed data points, as 
in figure 2. 
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their likelihoods are 0.034, 1.9 -10~ 3 , 0.45 -10~ 5 and 2.2 
•10~ 16 respectively. If we overestimate the rate of intrinsi- 
cally strong GRBs by a factor 2.5 (assuming that 4 of the 
observed events with Pj > 40 are a fluctuation, so that 
ho = 4, then the likelihoods are 0.40, 0.015, 0.95-10" 3 , 
2.65-10 -9 respectively. 

Note that the rejection of the NE model (10 -16 ) is now 
much stronger than in STS (10~ 4 ). This improvement is 
partially due to the better statistics of the joint BATSE- 
Ulysses sample, but mainly due to the fact that STS used 
too low an estimate for /40: 3 events per year, inferred 
from the conservative assumption that F s ~ 0.5 (com- 
pared to the present estimate F s ~ 0.1 obtained using 
Ulysses data). 

The sharp break in figure 5 around P a ~ 30 photons s _1 
cm~ 2 is statistical in origin (note the corresponding devi- 
ation of the two data points in figure 4). This fluctuation 
is mainly due to the Ulysses sample; it is evident in Atteia 
et al. (1999). 
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Fig. 5. — The same curves as in figure 4, but in integral form. 
Histogram - the integral peak flux distribution of Ulysses /BATSE 
GRBs. Solid curve: - SF model; dotted curve: NE; dashed curve: 
L95; long dashed curve: PZY98. 

Adopting a rate of intrinsically bright GRBs h = 10, 



Fig. 6. — The likelihood versus the fall-off factor (a + 1) where 
a is the parameter in equqtion 2. The arrow shows the result for 
the SF model; the 90% confidence limit (dashed horizontal line) is 
given with respect to this model. 

Figure 6 shows the likelihood factor for the parametrized 
source evolution model (equation 2) versus the fall-off fac- 
tor a+1. The results for NS merger models are also shown; 
the ordinate for these models is just the ratio of the max- 
imal NS merging rate (at z ~ 2) to that at z = 0. The 
likelihood curve has no turover at large a because our lu- 
minosity fuction has only a lower limit constraint at its 
bright end. It is interesting that the curve still displays a 
considerable increase (by a factor of 2.6) from a=22, which 
corresponds to the SF curve, to a = 80, i.e., the data are 
better fit by a GRB progenitor fall-off which is faster than 
the SF rate. This could be a natural consequence if the 
progenitors are supermassive stars whose population can 
decline faster than the total SF. However, this indication is 
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statistically weak and could also result, for example, from 
the same fluctuation in the Ulysses data which produces 
a break in the log N - log P curve around 20 - 30 photons 
s _1 cm~ 2 , as discussed above. 

More details of the fits are presented in Table f , where 
X 2 and the predictions for A 50 are given separately. We do 
not present the best fit parameters for the broken power 
law luminosity function because they are consistent with 
the results of STS where this issue was studied in detail. 

Table f , as well as figures 4 and 5, clearly demonstrate 
that models with an insufficiently steep evolutionary de- 
cline of GRB progenitors predict an excess of apparently 
strong GRBs with respect to the observations. 



Model 


lkh 


x 2 


^50 


lkh(A 50 ) 


NE 


2.2-f0- ib 


83 


49 


1.28-10 -a 


SF22 


0.034 


31 


18.6 


0.24 


L85 


1.9-10-5 


37 


33 


2.M0~ 4 


PZY98 


1.9-10" 3 


32 


26 


1.4-10" 2 


SF 80 


0.088 


31 


15 


0.57 


SF40 


0.062 


31 


16 


0.47 


SF 10 


0.36-10~ 2 


31 


25 


2.2-10~ 2 


SF 5 


0.76-10- 4 


38 


30 


1.9-10" 3 



Table 1 The maximum likelihood results for various mod- 
els. The second column (lkh) gives the final likelihood factor; 
the third column, the x 2 value (at 24 degrees of freedom); the 
fourth, the predicted A50 (the observed A 5 q is 15); the fifth, 
the probability of observing A50 less than 16 for its predicted 
value. The subscripts in the first column correspond to the 
value of a (see equation 2). SF22 in row 2 corresponds to the 
measured star formation curve. 

If we use the Bayesian approach, treating the ratio of 
likelihoods as the relative probabilities of different models, 
then the rejection factors for NS models relative to SF are 
0.055 for PZY98 and 1.3-10" 3 for L95. If we adopt the es- 
timate I40 — 4 then the constraints relax to 0.37 and 0.024 
respectively, i.e. the PZY98 model is consistent with the 
data. Note however that the choice of I40 = 4 corresponds 
to a less than 0.1 probablity fluctuation in the number of 
intrinsically bright GRBs with measured redshifts. 

The minimal fall-off factor allowed by the data at 90% 



confidence level is a + 1 ~ 12 or a + 1 ~ 8 for J40 = 4. 

5. CONCLUSIONS 

The joint BATSE - Ulysses data confirm a sharp de- 
cline in the GRB source population between z ~ 2 and 
the present epoch. Although it is consistent with that of 
star formation, a faster decline is slightly preferable, al- 
beit at a statistically insignificant level (~ la). The two 
models of binary system evolution leading to a final NS 
merger are well beyond the 90% confidence limit, except 
for the /40 = 4 case, which is based on the assumption 
of a large fluctuation in the observed number of intrinsi- 
cally bright GRBs. Note that while the statistics of bright 
GRBs will improve slowly, the redshift statistics can im- 
prove much faster, so that a more reliable estimate of I40 
may be available relatively soon. 

The joint BATSE/ Ulysses data present a new challenge 
to the neuton star binary model as an explanation of the 
source of long GRBs. Together with the results of after- 
glow studies it makes it very improbable. The only way 
to save the NS model is to show that the typical lifetimes 
of such systems is short. If very few survive longer than 
1 Gyr, this will fit the log N - log P distribution, and if 
many merge in a few Myr, this will explain the locations of 
the observed afterglows in the star forming regions of their 
host galaxies. Such a possibility has been studied in the 
recent work of Belczynski, Bulik & Rudak (2001) where 
it is shown that this could occur in some binary evolution 
models due to common envelope events producing very 
tight NS systems. Finally it should be pointed out again 
that our constraints refer only to the class of long GRBs, 
while the NS binary model is probably able to explain the 
origin of short bursts. 
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